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Abstract. We give a new probabilistic algorithm for interpolating a 
"sparse" polynomial / given by a straight-line program. Our algorithm 
constructs an approximation /' of /, such that / — /* probably has 
r )' at most half the number of terms of /, then recurses on the difference 

f — f ■ Our approach builds on previous work by Garg and Schost (2009), 

and Gicsbrecht and Roche (2011), and is asymptotically more efficient 

Y\ ' in terms of the total cost of the probes required than previous methods, 

in many cases. 



I> ! 1 Introduction 



We consider the problem of interpolating a sparse, univariate polynomial 

/ = ciz"^ + C2Z^^ +■■■ + ztz"^ £ n[z] 



^. .... 

f^ _ of degree d with t non-zero coefhcients ci , . . . , ct (where t is called the sparsity 

m ■ of /) over a ring TZ. More formally, we are given a straight-line program that 

evaluates / at any point, as well as bounds D > d and T > t. The straight- 
line program is a simple but useful abstraction of a computer program without 

^,j _ branches, but our interpolation algorithm will work in more common settings of 

'Oj ■ "black box" sampling of /. 

C^ [ We summarize our final result as follows. 

Theorem 1 Let f £ 7?.[z], where TZ is any ring. Given any straight-line program 
of length L that computes j , and bounds T and D for the sparsity and degree 
of f, one can find all coefficients and exponents of f using 0~{LT log D -\- 
LTlog D log(l//Lt))T ring operations in TZ, plus a similar number of bit operations. 
The algorithm is probabilistic of the Monte Carlo type: it can generate random 
bits at unit cost and on any invocation returns the correct answer with probability 
greater than 1 — /i, for a user supplied tolerance /i > 0. 



' For summary convenience we use soft-Oh notation: for functions (j>,ip £ R>o 
we say 4> € 0~{ip) if and only ii <j) £ 0{ip{logipy) for some constant c > 0. 



1.1 The straight-line program model and interpolation 

Straight-line programs are a useful model of computation, both as a theoretical 
construct and from a more practical point of view; see, e.g., (Biirgisser and Michael Clausen, 
1997, Chapter 4). Our interpolation algorithms work more generally for N- 
variate sparse polynomials / S 7?,[zi, . . . , zn] given by a straight-line program 
<S/ defined as follows. Sf takes an input {ai,. . . ^qn) G Ti-^ of length N, and 
produces a vector b G 7?.^ via a series of L instructions Ti : 1 < i < L oi the 
form 

_ j-fi < — ai-ka2, or 

I 7i -i — S E TZ (i.e., a constant from TZ), 

where • is a ring operation '-|-', ' — ', or 'x', and either ai,a2 G {ak}i<k<n 
or ai,a2 G {lj}i<j<i- When we say Sf computes /, we mean Sf sets 7l to 
/(ai,...,aAr) G Tl. 

To interpolate an A^- variate polynomial / G TZ[zi, . . . , zn], we apply a Kro- 
necker substitution, and interpolate 

m - / {z, z(^+^\ ^(^+^)^ . . . , z(^+i)"-^) G n[z]. 

While this certainly increases the degree, / and / have the same number of non- 
zero terms, and / can be easily recovered from /. This reduces the problem of 
interpolating the A'^- variate polynomial / of partial degree at most D to interpo- 
lating a univariate polynomial / of degree at most {D + 1)^ . For the remainder 
of this paper we thus assume / is univariate. 

It will also be necessary to evaluate our polynomial / G Tl[z], or rather 
our straight-line program Sf for /, in an extension ring of TZ. Precisely, we 
want to evaluate / at symbolic fth roots of unity for various choices of £, or 
algebraically, in 7?.[z]/(z^— 1). This may be regarded as transforming our straight- 
line program by substituting operations in TZ with operations in TZ[z]/{z — 1), 
where each element is represented by a polynomial in TZ[z] of degree less than 
i. Each instruction Ti in the transformed branching program now potentially 
requires M(^) operations in TZ, where M{£) is the number of operations in TZ and 
bit operations needed to multiply two degree-^ polynomials over the base ring 
TZ. By Cantor and Kaltofen (1991), we may assume M(^) = ^(^log^loglogf). 

Each evaluation of our straight-line program for / in TZ[z]/{z^ — 1) is called a 
probe of degree i. Thus, the cost of a degree-^ probe to Sf is 0~{L£) operations 
in TZ, and similarly many bit operations. 

This is easily connected to the more "classical" view of sparse interpolation, 
in which probes are simply evaluations of a "black box" polynomial at a single 
point (and we do not have any representation for how / is calculated). Each probe 
in the straight-line program model can be thought of as evaluating / at all £th 
roots of unity in the classical model. Since we charge M(^) = 0~(i) operations in 
TZ for a degree £ probe in the straight-line program model, i.e., about £ times as 
much as a single probe, this is consistent with the costs in a classical model. We 
note that algorithms for sparse interpolation presented below could be stated in 



this classical model, though we find the straight-line program model convenient 
and will continue with it throughout this paper. 

1.2 Previous work 

Straight-line programs, or equivalently algebraic circuits or DAGs, are impor- 
tant both as a computational model and as a data structure for polynomial 
computation. Their rich history includes both algorithmic advances and practi- 
cal implementations (Kaltofcn, 1989; Sturtivant and Zhang, 1990; Bruno et al., 
2002). 

One can naively interpolate an / G Tl[z] given by a straight-line program us- 
ing a dense method, with D probes of degree 1. Prony's (1795) interpolation algo- 
rithm — see (Bcn-Or and Tiwari, 1988; Kaltofcn et al., 1990; Gicsbrccht ct al., 
2009) — is a sparse interpolation method that uses evaluations at only 2T pow- 
ers of a root of unity whose order is greater than D. However, in the straight-line 
program model for a general ring, this would require evaluating at a symbolic 
Dth root of unity, which would use at least f2{D) ring operations and defeat the 
benefit of sparsity. Problems with Prony's algorithm are also seen in the classical 
model in that the underlying base ring TZ must also support an efficient discrete 
logarithm algorithm on entries of high multiplicative order (which, for example, 
is not feasible over large finite fields). 

We mention two algorithms specifically intended for straight-line programs. 

The Garg-Schost deterministic algorithm. Garg and Schost (2009) de- 
scribe a novel deterministic algorithm for interpolating a multivariate polyno- 
mial / given by a straight-line program. Their algorithm entails constructing an 
integer symmetric polynomial with roots at the exponents of /: 



^\[{v-e^)^nvi 



which is then factored to obtain the exponents e^. 

Their algorithm first finds a good prime: a prime p for which the terms of / 
remain distinct when reduced modulo z^—1. We call such an image / mod (z^— 1) 
a good image. Such an image gives us the values e^ mod p and hence xiv) mod p. 

Example 2 For f = z^^ + z^ , 5 is not a good prime because f mod (z^ — 1) = 
2z^. We say z^'^ and z^ collide modulo z^ — 1. 7 is a good prime, as the image 
f{z) mod {z'^ — 1) = z^ + z^ has as many terms as f{z) does. 

In order to guarantee that we have a good prime, the algorithm requires 
that we construct the images / mod {zP — 1) for the first N primes, where N 
is roughly 0~{T^ log D). A good prime will be a prime p for which the image 
/ mod {z^ — 1) has maximally many terms, which will be exactly t. Once we 
know we have a good image we can discard the images / mod (z'' — 1) for bad 
primes g, i.e. images with fewer than t terms. We use the remaining images to 



construct xiu) ~ Y[i=iiy ~ ^0 ^ "^[u] by way of Chinese remaindering on the 
images xiv) niodp. 

We factor xiv) to obtain the exponents e^, after which we directly obtain the 
corresponding coefficients Ci directly from a good image. 

The algorithm of Garg and Schost (2009) can be made faster, albeit Monte 
Carlo, using the following fact by Giesbrecht and Roche (2011). 

Fact 3 (Giesbrecht and Roche, 2011) Let f E 7?.[z] be a polynomial with at 
most T terms and degree at most D. Let X = niax(21, \^T(T — l)log_D]). A 
prime p chosen at random in the range [A, 2A] is a good prime for f{z) with 
probability at least i. 

Thus, in order to find a good image with probability at least 1 — e, we can 
inspect images / mod (z^ — 1) for [log 1/e] primes p chosen at random in [A, 2A]. 
As the height of xiv) can be roughly as large as D^ , we still require some 
C{T\ogD) probes to construct xiv)- 



The "diversified" interpolation algorithm. Giesbrecht and Roche (2011) 
obtain better performance by way of diversification. A polynomial / is said to 
be diverse if its coefficients Ci are pairwise distinct. The authors show that, for 
/ over a finite field or C and for appropriate random choices of a, f{az) is 
diverse with probability at least i. They then try to interpolate the diversified 
polynomial f{az). 

Once we have t with high probability, we look at images f{az) mod (z^—l) for 
primes p in [A, 2A], discarding bad images. As f{az) is diverse, we can recognize 
which terms in different good images are images of the same term. Thus, as 
all the Ci are at most £*, we can get all the exponents Ci by looking at some 
©"(log-D) good images of /. 



1.3 Deterministic zero testing 

Both the Monte Carlo algorithms of Garg and Schost (2009) and Giesbrecht and Roche 
(2011) can be made Las Vegas (i.e., no possibility of erroneous output) by way of 
deterministic zero-testing. Given a polynomial / represented by a straight-line 
program, each of these algorithms produces a polynomial /* that is probably /. 

Fact 4 (Blaser et al. (2009); Lemma 13) LetTZ be an integral domain. Sup- 
pose f — f* mod {zP — 1) for TlogD primes. Then f = f* ■ 

Thus, testing the correctness of the output of a Monte Carlo algorithm re- 
quires some C"'(T log D) probes of degree at most 0~{T log D). This cost does 
not dominate the cost of either Monte Carlo algorithm. We note that this deter- 
ministic zero test can dominate the cost of the interpolation algorithm presented 
in this paper if T is asymptotically dominated by logl^. 



1.4 Summary of results 

We state as a theorem the number and degree of probes required by our new 
algorithm presented in this paper. 

Theorem 5 Let f G TZ[z], where TZ is a ring. Given a straight-line program for 
f, one can find all coefficients and exponents of f with probability at least I — jJ, 
using O" ( logT(logiD + log i) 1 probes of degree at most 0(T \og^ D). 



Table 1. A "soft-Oh" comparison of interpolation algorithms for straight-line programs 





Probes 


Probe degree 


Cost of probes 


Type 


Dense 


D 


1 


LD 


deterministic 


Garg & Schost 


T'-'logD 


T^logD 


LT-" log^ D 


deterministic 


*Las Vegas G & S 


TlogD 


T' log D 


LT-" log" D 


Las Vegas 


* Diversified 


\ogD 


T' log D 


LT^ log" D 


Las Vegas 


fRecursive 


log TlogD 


Tlog^D 


LT log-" D 


Monte Carlo 



^Average # of probes given; f for a fixed probability of failure /^ 



Table 1 gives a rough comparison of known algorithms. Our recursive al- 
gorithm improves by a "soft-Oh" factor of T/logD over the Giesbrecht- Roche 
diversification algorithm, and as such is better suited for moderate values of 
T. Our algorithm recursively interpolates a series of polynomials of decreasing 
sparsity. An advantage of this method is that, when we cross a threshold where 
logD begins to dominate T, we can merely call the Monte Carlo diversification 
algorithm instead. 



2 A recursive algorithm for interpolating / 



Entering each recursive step in our algorithm we have our polynomial / repre- 
sented by a straight-line program, and an explicit sparse polynomial /* "approx- 
imating" /, that is, whose terms are generally a subset of the terms in the sparse 
representation of /. At each recursive step we try to interpolate the difference 
g = f — f* ■ To begin with, /* is initialized to zero. 

We first find an "ok" prime p which separates most of the terms of g. We 
then use that prime p to build a approximation /**, containing most of the 
terms of g, plus possibly some additional "deceptive" terms. The polynomial 
/** is constructed such that g ~ f ~ f* has, with high probability, at most T/2 
terms. We then recursively interpolate the difference g — f** ■ 

Producing images /* mod {z^ — 1) is straightforward, we merely reduce the 
exponents of terms of /* modulo L We assume g has a sparsity bound Tg < T. 



2.1 A weaker notion of "good" primes 

To interpolate a polynomial g, the sparse interpolation algorithm described by 
Giesbrecht and Roche (2011) requires a good prime p which keeps the exponents 
of g distinct modulo p. That is, g mod (z^ — 1) has the same number of terms 
as g. We define a weaker notion of a good prime, an ok prime, which separates 
most of the terms of g. To that end we measure, for fixed g and prime p, how 
well p separates the terms of g. 

Definition 6 Fix a polynomial g = J2i=i '^i^^^ ^ ^[^] ''^'''^^ non-zero ci, . . . ,Ct G 
TZ, where Ci < Cj for i < j, we say Ciz'^' and CjZ'^^ , i ^ j, collide modulo z^ — 1 
if Bi = Cj mod p. We call any CiZ*^^ of f which collides with any other a colliding 
term of f modulo z^ — 1. We let Cg{p) G [0,t] denote the number of colliding 
terms of g modulo z'p — 1. 

Example 7 For the polynomial g — 1 + z'"' + 2:^ + z^", Cg(2) = 4, since 1 collides 
with z^^ and z° collides with z'^ modulo z^ — 1. Similarly, Cg(5) — 2, since z^ 
collides with z^'^ modulo z^ — 1. 

We say CiZ^^ and Cjz"^^ collide modulo z^ — 1 because both terms have the 
same exponent once reduced modulo z^ — 1. All other terms of g we will call 
non- colliding terms modulo z^ — 1. 

In the sparse interpolation algorithm of Giesbrecht and Roche (2011), one 
chooses a A G Z>o such that the probability of a prime p G [A, 2A], chosen at 
random, having Cg{p) = 0, is at least i However, in order to guarantee that we 
find such a prime with high probability, we need to choose A G 0{T^ logD). 

In this paper we will search over a range of smaller primes, while allowing 
for a reasonable number of collisions. We try to pick A such that 

Pr {Cg{p) > 7 for a random prime p G [A, 2A]) < 1/2, 

for a parameter 7 to be determined. 

Lemma 8 Let g £ TZ[z] be a polynomial with t < T terms and degree at most 

d < D. Suppose we are given T and D, and let A = max I 21, ^ ~ -^ "^ — - ] . 

Let p be a prime chosen at random in the range A, . . . , 2A. Then Cg{p) > 7 with 
probability less than i. 

Proof. The proof follows similarly to the proof of Lemma 2.1 in (Giesbrecht and Roche, 
2011). 

Let B be the set of unfavourable primes for which Cg{p) > 7 terms collide 
modulo zP — 1, and denote the size of B by #_B. As every colliding term collides 
with at least one other term modulo z^ — 1, we know p^sip) divides ni<i=^7<t(s«~ 
Cj). Thus, as Cg{p) > 7 for p G -B, 

X*^! < J] p^ < [] (e, - e,) < d*(*-i) < 7?^(^-i). 



Solving the inequality for #5 gives us 

T{T~ l)ln(D) 



#B< 



ln(A)7 



The total number of primes in [A, 2A] is greater than 3A/(51n(A)) for A > 21 by 

Corollary 3 to Theorem 2 of (Rosscr and Schoenfeld, 1962). From our definition 

of A we have 

3A ^ 2TiT-l)HD) 

51n(A) ln(A)7 "^ ' 

completing the proof. D 

Relating the sparsity of g mod (z^ — 1) virith Cg{p) 

Suppose we choose A according to Lemma 8, and make k probes to compute 
g mod {zP^ — 1), . . . , g mod (z^* — 1). One of the primes pi will yield an image 
with fewer than 7 colliding terms (i.e. Cg{pi) < 7) with probability at least 
1 — 2~^. Unfortunately, we do not know which prime p maximizes Cg{p). A good 
heuristic might be to select the prime p for which g mod (z'' — 1) has maximally 
many terms. However, this does not necessarily minimize Cg{p). For instance, 
consider the following example. 



Example 9 Let 

We have 



l + z + z^-2z 



4 0^13 



g mod (z — 1) = 2 — z, and g mod (z — 1) = 1. 

While g mod (z^ — 1) has more terms than g mod {z^ — 1), we see that Cg{2) — 4 
is larger than Cg{3) ~ 3. 

While we cannot determine the prime p for which g mod (z^ — 1) has maxi- 
mally many non-colliding terms, we show that choosing the prime p which max- 
imizes the number of terms in g mod (z^ — 1) is, in fact, a reasonable strategy. 

We would like to find a precise relationship between Cg{p), the number of 
terms of g that collide in the image g mod (z^ — 1), and the sparsity s of 5 mod 
(zP-1). 

Lemma 10 Suppose that g has t terms, and g mod (z^ — 1) has s < t terms. 
Then t- s < Cg{p) < 2{t ~ s). 

Proof. To prove the lower bound, note that t ~~ Cg {p) terms of g will not collide 
modulo z^ — 1, and so g mod (z^ — 1) has sparsity s at least t — Cg{p). 

We now prove the upper bound. Towards a contradiction, suppose that 
Cg{p) > 2{t — s). There are Cg{p) terms of g that collide modulo z^ — 1. Let 
h be the Cg (p)-sparse polynomial comprised of those terms of g. As each term of 
h collides with at least one other term of h, h mod (z^ — 1) has sparsity at most 
Cg{p)/2. Since none of the terms oig — h collide modulo z^ — 1, g — h mod (z^ — 1) 



has sparsity exactly t — Cg{p). It follows that g mod (z^ — 1) has sparsity at most 
t-Cg{p)+Cg{p)/2 = t-Cg{p)/2. That is, s < t-Cg{p)/2, andsoCg(p) < 2{t~s). 

a 

Corollary 11 Suppose g mod (z* — 1) has sparsity s^, and g m.od (z^ — 1) has 
sparsity Sp > Sg. Then Cg{p) < 2Cg{q). 

2 roof 

Cg{p) < 2{t — Sp) by the second inequality of Lemma 10, 

< 2(t — Sq) since Sp > Sq, 

< 2Cg{q) by the first inequality of Lemma 10. D 

Suppose then that we have computed g mod {z^ — 1), for p belongmg to some 
set of primes S, and the minimum value of Cg{p), p € S, is less than 7. Then 
a prime p* G S* for which g mod (z^ — 1) has maximally many terms satisfies 
Cg{p*) < 27. We will call such a prime p* an ok prime. 

We then choose 7 = wT for an appropriate proportion w e (0, 1). We show 
that setting w = 3/16 allows that each recursive call reduces the sparsity of 
the subsequent polynomial by at least half. This would make A — [— 10/3(T — 
l)ln(r')] = \^{T — l)ln(Z3)]. As per Lemma 8, in order to guarantee with 
probability 1 — e that we have come across a prime p such that Cg{p) < 7, we will 
need to perform [log 1/e] probes of degree 0{TlogD). Procedure FindOkPrime 
reiterates how we find an ok prime. 



Procedure FindOkPrime (Sf, f*,Tg,D, e) 

Input: 

— Sf, & straight-line program that computes a polynomial / 

— /*, a current approximation to / 

— Tg and D, bounds on the sparsity and degree oi g = f — f* respectively 

— £, a bound on the probability of failure 

Output: With probability at least 1 — e, we return an "ok prime" for <? = / — /* 

X< — max (21, \if-{Tg-l)lnD]) 
max_sparsity,p i — 0,0 
for i < — 1 to [log 1/e] do 

p' < — a random prime in [A, 2A] 

if # of terms of (/ — /*) mod [z'' — 1) > max.sparsity then 
max_sparsity < — # of terms of (/ — /*) mod {z^ — 1) 
pi — p' 

return p 



A practical application would probably choose random primes by selecting 
random integer values in [A, 2A] and then applying probabilistic primality test- 
ing. In order to ensure deterministic worst-case run-time, we could pick random 
primes in the range [A, 2A] by using a sieve method to pre-compute all the primes 
up to 2A. 



2.2 Generating an approximation /** of g 

We suppose now that we have, with probabihty at least 1 — e, an ok prime p, 
i.e., a prime p such that Cg{p) < 2wT for a suitable proportion w. We now use 
this ok prime p to construct a polynomial /** containing most of the terms of 

g = .!-!*■ 

For a set of coprime moduli Q = {gi, . . . , q^} satisfying nj=i 9* > ^i ^^ will 
compute g mod (z^''' — 1) for q < i < k. Here we make no requirement that the 
qi be prime. We merely require that the qi are pairwise co-prime. 

We choose the qi as follows: denoting the i"^ prime by pi, we set qi = Pi ^^ , 
for an appropriate choice of x. That is, we let qi be the greatest power of the i"* 
prime that is no more than x. For pi < x, we have qi > x/pi and qi > pi Either 
x/pi or Pi is at least ^/x, and so qi > \fx as well. 

By Corollary 1 of Theorem 2 in Rosscr and Schoenfeld (1962), there are more 
than a;/ In a; primes less than or equal to x for x > 17. Therefore 



n * > (v^) 

Pi<x 



X I In a 



As we want this product to exceed _D, it suffices that 
lui?<ln((Vi)"^'"^)=x/2. 

Thus, if we choose x > max(21n(Z?), 17), and k = [a;/lna;], then ni=i 9j ^^^^ 
exceed D. This means qi € 0{\ogD) and pqi G 0{Tlog^ D). The number of 
probes in this step is fc € 0(log(Z?)/loglog(I?)). Since we will use the same set 
of moduli Q = {gi, . . . , qk} in every recursive call, we can pre-compute Q prior 
to the first recursive call. 

We now describe how to use the images g mod (z^'?' — 1) to construct a 
polynomial /** such that g — /** is at most r/2-sparse. 

If cz^ is a term of g that does not collide with any other terms modulo z^ — 1, 
then it certainly will not collide with other terms modulo z^* ~ 1 for any natural 
number q. Similarly, if c*z'^ modp appears in g mod {z^ — 1) and there exists a 
unique term c* z'^ ™°'* p^' for i = 1, 2, . . . k, then c*z'^ is potentially a term of g. 
Note that c*z'^ is not necessarily a term of g: consider the following example. 

Example 12 

giz) ^l + z + z^ + z^ + z"+4 - zi4-ii+4 _ ^15.11+4^ 

with hard sparsity bound Tg ~ 7 and degree bound D = 170 and let p — 11. We 
have 

g{z) mod (z" - I) = I + z + z"^ + z^ ~ z^. 



10 



As dcg(g) = 170 < 2 • 3 • 5 • 7 = 210, it suffices to make the probes g mod z^^'' — 1 
for q = 2,3, 5, 7. Probing our remainder black-box polynomial, we have 



g mod (z^^ - 1) = 1 
g mod (z^^ - 1) = 1 
g mod (z^^ - 1 
g m.od {z — 1 



z3-zl5 



z3 - z26, 



1 + z + z^ + z^- z^\ 
1 + z + z^ + z^ - z^'l 



/n eac/i of the images g mod z'"' — 1, t/iere is a unique term whose degree is 
congruent to one of e ^ 0, 1, 2, 3, 4 modulo p. Four of these terms correspond to 
the terms l,z,z^,z*^ appearing in g. Whereas the remaining term has degree e 
satisfying e = 1 mod 2, e = 2 mod 3, e = 3 mod 5, and e = 1 mod 7. By Chinese 
remaindering on the exponents, this gives a term —z^^^ not appearing in g. 



Procedure ConstructApproximation(5/, /*,Tg, D,p, Q) 

Input: 

— Sf, & straight-line program that computes a polynomial / 

— /*, a current approximation to / 

— Tg and D, bounds on the sparsity and degree oi g = f — f* respectively 

— p, an ok prime for g (with high probabihty) 

— Q, a set of co-prime moduli whose product exceeds D 

Output: A polynomial /** such that, if p is an ok prime, g — /" has sparsity 
at most [Tg/2]. 



II Collect images of g 

£ < — set of exponents of terms in (f — f*) mod (z^ 

for q (z Q do 

h^{f-r)mod{z^^~l) 

for each term cz'^ in h do 



if -B(emodp).ij is already initialized then £ < — £/{e mod p} else 
E(e mod p),q < — « mod q 



II Construct terms of new approximation of g, /** 

for Cp £ £ do 

e i — least nonnegative solution to {e — Eej,,q mod g | g G Q} 
c i — coefficient of z^'' term in (f — f*) mod (z"^ — 1) 
if e < D then /" i — /" + cz"' 

return /** 



Definition 13 Let c*z^ be a monomial such that c*z'^ modp appears in g m.od 
zP — 1, and c*z'^ mod gi jg ^^g unique term of degree congruent to e* modulo p 



11 

appearing in g mod (z^^' — 1) for each modulus qi. If c*z'^ is not a term of g we 
call it a deceptive term. 

Fortunately, we can detect a collision comprised of only two terms. Namely, 
if ciz'^'^ + C2z'^^ collide, there will exist a qi such that qi | (ei — 62). That is, 
g mod (zP''^ — 1) will have two terms whose degree is congruent to ei mod p. Once 
we observe that, we know the term (ci +C2)z'^'^ "^"'^ ^ appearing in g mod (z^ — 1) 
was not a distinct term, and we can ignore exponents of the congruence class 
ei mod p in subsequent images g mod (z^'^^ — 1). 

Thus, supposing g mod {z^ — 1) has at most 2j colliding terms and at least 
t — 27 non-colliding terms, /** will have the i — 27 non-colliding terms of g, plus 
potentially an additional |7 deceptive terms produced by the colliding terms of 
g. In any case, g — f** has sparsity at most I7. Choosing 7 = j^t guarantees that 
g-f** has sparsity at most t/2 < T/2. This would make A = \^{T- 1) \n{D)'] . 

Procedure Construct Approximation gives a pseudocode description of how 
we construct /**. 

If we find a prospective term in our new approximation /** has degree greater 
than D, then we know that term must have been a deceptive term an discard it. 
There are other obvious things we can do to recognize deceptive terms which we 
exclude here. For instance, we should check that all terms from images modulo 
zP'i — 1 whose degree agree modulo p share the same coefficient. 



2.3 Recursively interpolating f — f* 



Procedure Interpolate (5/, T, D,£) 

Input: 

— Sf, a, straight-line program that computes a polynomial / 

— T and D, bounds on the sparsity and degree of /, respectively 

— fj,, an upper bound on the probability of failure 
Output: With probability at least 1 — ^, we return / 

xi — max(21n(D),17) 

Q i — {p^^°^p ^J : p is prime, p < x} 

return InterpolateRecurse(5/,0, T, D, Q,^/(logr-|- 1)) 



Once we have constructed /**, we refine our approximation /* by adding /** 
to it, giving us a new difference g = / — /* containing at most half the terms 
of the previous polynomial g. We recursively interpolate our new polynomial g. 
With an updated sparsity bound lTg/2\ , we update the values of 7 and A and 
perform the steps of Sections 2.1 and 2.2. We continue in this fashion for logT 
times. Thus, the total number of probes becomes 



o(\ogT{j 



logo 
>g log D 



log(l/e)) 
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Procedure InterpolateRecursedS/, /*,Tg,D, Q,e) 

Input: 

— Sf, & straight-line program that computes a polynomial / 

— /*, a current approximation to / 

— Tg and D, bounds on the sparsity and degree oi g = f ~ f* , respectively 

— Q, a set of coprime moduli whose product is at least D 

— e, a. bound on the probability of failure at one recursive step 
Output: With probability at least 1 — fJ,, the algorithm outputs / 

if Tg — then return /* 

p< — FindOkPrime(Sf, f ,Tg,V,e) 

/** i — ConstructApproximation{Sf,f*,Tg,D,p,Q) 

return InterpolateRecurse(5/, /* +/**, \Tg\,D, Q,e) 



of degree at most 0(riog^ D). 

Note now that in order for this method to work we need that, at every 
recursive call, we in fact get a good prime, otherwise our sparsity bound on 
the subsequent difference of polynomials could be incorrect. At every stage we 
succeed with probability 1 — £, thus the probability of failure is 1 — (1 — e) ri°sri _ 
This is less than [logr]£. If we want to succeed with probability fi, then we can 
choose e = i^;^eO(i^). 

Interpolate pre-computes our set of moduli Q, then makes the first recur- 
sive call to InterpolateRecurse, which subsequently calls itself. 

2.4 A cost analysis 

We analyse the cost of our algorithm, thereby proving Theorems 1 and 5. 

Pre-computation. Using the wheel sieve (Pritchard, 1982), we can compute 
the set of primes up to a; S 0{\ogD) in 0~{logD) bit operations. From this set 
of primes we obtain Q by computing pL'ogp^^J for p < ,y^ by way of squaring- 
and-multiplying. For each such prime, this costs O~{\ogx) bit operations, so the 
total cost of computing Q is 0~{\ogD). 

Finding ok primes. In one recursive call, we will look at some logl/e = 
©(log 1/^ log log T) primes in the range [A, 2A] in order to find an ok prime. 
Any practical implementation would select such primes by using probabilistic 
primality testing on random integer values in the range [A, 2A]; however, the 
probabilistic analysis of such an approach, in the context of our interpolation al- 
gorithm, becomes somewhat ungainly. We merely note here that we could instead 
pre-compute primes up to our initial value of A G 0{TlogD) in 0~{T log D) bit 
operations by way of the wheel sieve. 

Each prime p is of order T log D, and so, per our discussion in Section 1, each 
probe costs 0~{LT log D) ring operations and similarly many bit operations. 
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Considering the 0{\ogT) recursive calls, this totals C'"'(LT log L* log l//i) ring 
and bit operations. 



Constructing the new approximation /**. Constructing /** entails ©"(log-D) 
probes of degree 0~{T \og^ D). This costs 0"'(LTlog^ D) ring and bit operations. 
Performing these probes at each ©(logT) recursive call introduces an additional 
factor of log T, which does not affect the "soft-Oh" complexity. This step domi- 
nates the cost of the algorithm. 

Building a term cz"^ of /**, by Theorem 5.8 of Gathcn and Gerhard (2003), 
requires some O(log D) word operations. Thus the total cost of Chinese remain- 
dering to construct /** becomes 0{T log D). Again, the additional logT factor 
due to the recursive calls does not affect the stated complexity. 



3 Conclusions and Future Work 

We have presented a recursive algorithm for interpolating a polynomial / given 
by a straight-line program, using probes of smaller degree than in previously 
known methods. We achieve this by looking for "ok" primes which separate 
most of the terms of /, as opposed to "good" primes which separate all of the 
terms of /. As is seen in Table 1, our algorithm is an improvement over previous 
algorithms for moderate values of T. 

This work suggests a number of problems for future work. We believe our 
algorithms have the potential for good numerical stability, and could improve on 
Giesbrecht and Roche's (2011) work on numerical interpolation of sparse com- 
plex polynomials, hopefully capitalizing on the lower degree probes. Our Monte 
Carlo algorithms are now more efficient than the best known algorithms for 
polynomial identity testing, and hence these cannot be used to make them error 
free. We would ideally like to expedite polynomial identity testing of straight- 
line programs, the best known methods currently due to Blaser et al. (2009). 
Finally, we believe there is still room for improvement in sparse interpolation 
algorithms. The vector of exponents of/ comprises some TlogD bits. Assuming 
no collisions, a degree-^ probe gives us some tlogH bits of information about 
these exponents. One might hope, aside from some seemingly rare degenerate 
cases, that logD probes of degree Tlog-D should be sufficient to interpolate /. 
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